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Abstract 

In this paper, dynamical systems theory and bifurcation theory are applied to investi¬ 
gate the rich dynamical behaviours observed in three simple disease models. The 2- and 
3-dimensional models we investigate have arisen in previous investigations of epidemiol¬ 
ogy, in-host disease, and autoimmunity. These closely related models display interesting 
dynamical behaviors including bistability, recurrence, and regular oscillations, each of 
which has possible clinical or public health implications. In this contribution we elucidate 
the key role of backward bifurcation in the parameter regimes leading to the behaviors 
of interest. We demonstrate that backward bifurcation facilitates the appearance of Hopf 
bifurcations, and the varied dynamical behaviors are then determined by the properties of 
the Hopf bifurcation(s), including their location and direction. A Maple program devel¬ 
oped earlier is implemented to determine the stability of limit cycles bifurcating from the 
Hopf bifurcation. Numerical simulations are presented to illustrate phenomena of interest 
such as bistability, recurrence and oscillation. We also discuss the physical motivations 
for the models and the clinical implications of the resulting dynamics. 
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3 Introduction 


In the mathematical modelling of epidemic diseases, the fate of the disease can be predicted 
through the uninfected and infected equilibria and their stability. The basic reproduction num¬ 
ber, Rq, represents the average number of new infectives introduced into an otherwise disease- 
free system by a single infective, and is usually chosen as the bifurcation parameter. If the 
model involves a forward bifurcation, the uninfected equilibrium is in general globally asymp¬ 
totically stable Il28l . characterized hy Rq < 1, and infection fails to invade in this parameter 
regime. The threshold = 1 defines a bifurcation (or critical) point, and when Rq > 1, a 
stable infected equilibrium emerges. This simple exchange of stability implies that complex 
dynamics will not typically occur in forward bifurcation. 

In contrast, backward bifurcation describes a scenario in which a turning point of the in¬ 
fected equilibrium exists in a region where all state variables are positive, and < 1 • This 
induces multiple infected equilibria, disrupting the global stability of the uninfected equilib¬ 
rium. Multiple stable states (e.g., bistability) may likewise appear in [Il5l|4ll3, and Yu et al. 
(submitted for publication). Instead of converging globally to the uninfected equilibrium when 

< 1, the solution may approach an infected equilibrium, depending on initial conditions. 

In practice, the phenomenon of backward bifurcation gives rise to new challenges in dis¬ 
ease control, since reducing Rq such that i?o < 1 is not sufficient to eliminate the disease [l22l [5ll. 
Instead, Rq needs to be reduced past the critical value given by the turning point ll22ll . since 
the result in Yu et al. (submitted for publication) shows that the uninfected equilibrium in 
backward bifurcation is globally stable if Rq is smaller than the turning point. Furthermore, an 
infective outbreak or catastrophe may occur if Rq increases and crosses unity, while the upper 
branch of the infected equilibrium remains stable [fT5ll2Tll47ll4^ . In addition, oscillation or 
even recurrent phenomena may occur if uninfected and infected equilibria coexist in a param¬ 
eter range, and both are unstable [l47l [48l . [|23 predicted oscillations arising from backward 
bifurcation, and [HI pointed out that the unstable infected equilibrium “commonly arises from 
Hopf bifurcation”, but did not demonstrate oscillations. 

Several mechanisms leading to backward bifurcation have been proposed, such as partially 
effective vaccination programs [[5113, educational influence on infectives’ behavior [123 . the 
interaction among multi-group models flU [TOl and multiple stages of infection [|4^ . In 
this study, we will investigate the emergence of backward bifurcation in three simple disease 
models which have arisen in the study of epidemiology, in-host disease and autoimmunity. In 
each case, we find that backward bifurcation facilitates the emergence of Hopf bifurcation(s), 
and Hopf bifurcation in turn underlies a range of complex and clinically relevant dynamical 
behaviors. 

A central theme in our investigation is the role of the incidence rate in the epidemiological 
and in-host disease models. The incidence rate describes the speed at which an infection 
spreads; it denotes the rate at which susceptibles become infectives. Under the assumptions 
of mass action, incidence is written as the product of the infection force and the number of 
susceptibles. For example, if S and I denote the susceptible and infective population size 
respectively, a bilinear incidence rate, f{SJ) = SI (where /3 is a positive constant), is linear 
in each of the state variables: S and 7. 
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The possibility of saturation effects dH |71 has motivated the modification of the incidence 
rate from bilinear to nonlinear. Saturation occurs when the number of susceptible contacts 
per infective drops off as the proportion of infectives increases. A nonlinear incidence rate, 
therefore, typically increases sublinearly with respect to the growth of the infective population, 
and may finally reach an upper bound. The development of nonlinear incidence was first 
investigated in the form /3/^S^, where j5,p, and q are positive constants [[^l3ni2^l2^[T3ll29ll . 
Other forms of nonlinear incidence have also been analysed, such as kIPS/{ \ + al^) [|32l, and 
kS\n{l+vP/k) a. 

Since the nonlinear incidence functions described above were often developed to incor¬ 
porate saturation effects, these functions are typically concave at realistic parameter values. 
[1281 used this feature to derive general results for disease models with concave incidence. 
They proved that standard epidemiological models with concave incidence functions will have 
globally asymptotically stable uninfected and infected equilibria for < 1 and re¬ 

spectively. 

More specifically, denoting the incidence rate function as /(5, 7, N), where N is the popu¬ 
lation size, the classical SIRS model considered in [|2^ takes the form 

^=pN-f{SJ,N)-pS + aR, ^=/(5,/,A)-(5-fAt)/, ^ = dl-aR-pR, (1) 

at at at 

where jl, 5, and a represent the birth/death rate, the recovery rate and the loss of immunity 
rate, respectively. When a = 0, system ([T]) becomes an SIR model. Assuming that the total 
population size is constant, that is, N = S +1 + R, the above system can be reduced to a 2- 
dimensional model: 

«-l C T 

- = {a + p)N- /(S, l,N)-al-{a^p)S, - = f{S, l,N)-{5 + p)L (2) 

Moreover, it is assumed in [128ll that the function /(S, 7, N), denoting the incidence rate, satis¬ 
fies the following three conditions: 


/(S, 0,A)=/(0,7,A)=0, 


df{S, I,N) 
dl 

d^f{S, I,N) 
dl^ 


> 0 , 


dfiS,I,N) 


dS 

<0, V5,7>0. 


>0, V5,7>0 


(3a) 

(3b) 

(3c) 


The first two conditions (3a) and (3b) are necessary to ensure that the model is biologically 
meaningful. The third condition (3c) implies that the incidence rate f{S,I,N), is concave 
with respect to the number of infectives. It is also assumed that evaluated at the 

uninfected equilibrium is proportional to the basic reproduction number 7?o and thus 
should be a positive finite number [|2^ . Korobeinikov and Maini first considered 7 = 0, or 
/(5, 7, N) — {5 + fl)I = 0, and showed that forward bifurcation occurs in model ([^ with 
a concave incidence function. They further proved that the uninfected equilibrium Qq = 
(5o, 7o) = (A, 0) and the infected equilibrium Q= {S,I) are globally asymptotically stable. 


when Rq = 


1 df{So,lo,N) 
S+fl dl 


< 1 and7?o > 1, respectively. 
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In the sections to follow, for an incidence rate function /(5, 1), satisfying ( [^ and (3b), we 
define /(5, 1) as concave, if it satisfies (3c); as convex, if ^ > 0, V / > 0; and as convex- 

> 0, V/ e (0,/2), and > 0, 


di 


concave, if there exist 0 < I\ < h < +°o, such that 
V/e (0,/i), = 0, for 7 =/i, < 0, V/G (/i,/ 2 ). 

Several models closely related to ([^ have been previously studied. For example, by adding 


a saturating treatment term to model ([^ with a concave incidence rate, [|49l showed that this 
model may yield backward bifurcation and Hopf bifurcation. With an even more sophisticated 
nonlinear incidence rate function: kPS/{ \ + (xl^), where p = l = 2, [|3^ proved that a reduced 
2-dimensional SIRS model could exhibit backward bifurcation, Hopf bifurcation, and even 
Bogdanov-Takens bifurcation and homoclinic bifurcation. Although the choice of p = I = 
2 was not motivated by a specific physical process, this important result demonstrates that 
a nonlinear incidence rate can induce backward bifurcation, and further generate complex 
dynamics in a simple disease model. 


One of the focal points of our study will be a convex incidence function which arose in a 
4-dimensional HIV antioxidant therapy model ll43l . In this model, the infectivity of infected 
cells was proposed to be an increasing function of the density of reactive oxygen species, 
which themselves increase as the infection progresses. In P3l . meaningful parameter values 
were carefully chosen by data fitting to both experimental and clinical results. In this parameter 
regime, the model was observed to capture the phenomenon of viral blips, that is, long periods 
of undetectable viral load punctuated by brief episodes of high viral load. Viral blips have been 
observed clinically in HIV patients under highly active antiretroviral therapy [fTTl[T4ll35ll34ll . 
and have received much attention in the research literature, both by experimentalists lUTl [HI 
l20l and mathematicians [Ii6ii27i dal EH m. Nonetheless, the mechanisms underlying this 
phenomenon are still not thoroughly understood [|20l[3^ . 

We recently re-examined the model developed in ll43ll . with the aim of providing new 
insight into the mechanism of HIV viral blips [|47ll4^ . Focusing on the dynamics of the slow 
manifold of this model, we reduced the dimension of the 4-dimensional model by using quasi¬ 
steady state assumptions. After a further generalization and parameter rescaling process, a 
2-dimensional in-host HIV model [|47l l48ll was obtained, given by 


" = i 
dT 


AY 

- )XY, 


dY , AY 

— = (BH- 

dr ^ Y + C 


)XY-Y, 


(4) 


Y + C 

where X and Y denote the concentrations of the uninfected and infected cells respectively. The 
constant influx rate and the death rate of Y have been scaled to 1. The death rate of X is D. The 
2-dimensional infection model above (|^, reduced from the 4-dimensional HIV model [l43ll . 
preserves the viral blips observed in the HIV model. 

Importantly, system (|^ is equivalent to the SIR model (|^, except that the incidence func¬ 
tion is convex, as we will show in section [A2j This equivalence can be demonstrated if we set 


S = eix, I = e 2 y, and t = with ei = 62 = 


and €2 = 


rescaled to 


S+^' 


In this case, system (2) is 


dv ^ p 
dr d + p 




dy 1 ^ 


which takes the same form as system (Q. Therefore, although system ([^ arises in epidemi¬ 
ology and system (ffl was derived as an in-host model, they are mathematically equivalent in 
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this sense. We will refer to both systems Q and Q as infeetion models. 

In previous work P7ll48l . we analyze the recurrent behavior which emerges in system Q 
in some detail. Recurrence is a particular form of oscillatory behavior characterized by long 
periods of time close to the uninfected equilibrium, punctuated by brief episodes of high in¬ 
fection P31 . Thus HIV viral blips are an example of recurrent behavior, but recurrence is a 
more general feature of many diseases [l45ll4^ . We have demonstrated that the increasing and 
saturating infectivity function of system Q is critical to the emergence of recurrent behaviour. 
This form of an infectivity function corresponds to a convex incidence rate function in the as¬ 
sociated 2-dimensional infection model Q, and can likewise induce recurrence in this model. 
Convex incidence has been previously suggested to model ‘cooperation effects’ in epidemiol¬ 
ogy [|281 , or cooperative phenomena in reactions between enzyme and substrate, as proposed 
by [sa. 

The rest of this paper is organized as follows. In Section 2, we study two 2-dimensional 
infection models, both closely related to system Q. We show that system ([^ with either 
(a) a concave incidence rate and saturating treatment term or (b) a convex incidence rate as 
shown in system Q, can exhibit backward bifurcation; we then identify the necessary terms 
in the system equations which cause this phenomenon. In Section 3, we demonstrate that 
in both models, backward bifurcation increases the likelihood of a Hopf bifurcation on the 
upper branch of the infected equilibrium. Studying system Q in greater detail, we illustrate 
how the location of the Hopf bifurcations and their directions (supercritical or subcritical), 
determine the possible dynamical behaviors, concluding that backward bifurcation facilitates 
Hopf bifurcation(s), which then underly the rich behaviours observed in these models. In 
Section 4, we explore backward bifurcation further, presenting an autoimmune disease model 
which exhibits negative backward bifurcation, that is, a bifurcation for which the turning point 
when Rq < 1 is located in a region where one or more state variables is negative. Although 
this bifurcation introduces two branches of the infected equilibrium, we demonstrate that, in 
the biologically feasible area, only forward bifurcation exists in this model. We then present a 
modification to this autoimmune model, motivated by the recent discovery of a new cell type, 
which generates a negative backward bifurcation and Hopf bifurcation, and allows recurrent 
behavior to emerge. A conclusion is drawn in Section 5. 


4 Backward bifurcation 

In this section, we study backward bifurcation in two 2-dimensional infection models. In 
particular, we explore the essential terms and parameter relations which are needed to generate 
backward bifurcation. Furthermore, we examine the convex incidence rate, and reveal its 
underlying role in determining the emergence of backward bifurcation. 
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4.1 Backward bifurcation in the infection model with concave incidence 


First, we consider the SIR model with concave incidence, described by the following equa¬ 
tions 


d5 

dt 


A 


/ 35 / 

l+kl 


dS, 


dl 

dt 


psi 

1+kI 


{d Y 


dR 

dt 


yl - dR, 


( 5 ) 


where S, 1 and R denote the number of susceptible, infective, and recovered individuals, re¬ 
spectively; A is the constant recruitment rate of susceptibles; d, y, and £ represent the rates of 
natural death, recovery, and the disease-induced mortality, respectively. Note that the function 
is an incidence rate of the form lEHl, when I = h = 1. Here, /3 is the infection 

rate, and k measures the inhibition effect. Since the variable R is not involved in the first two 
equations, system Q can be reduced to a 2-dimensional model as 


dS 

dt 


A 


psi 

1+kI 


dS, 


dl 

dt 


psi 

l+kl 


{d + y+£)I. 


( 6 ) 


In [1^ , an additional assumption regarding limited medical treatment resources is introduced 
to the above model, leading to a model with a saturating treatment term, given by 


d5 

dt 


MS, I)=A 


psi 

l+kl 


— dS, 


dl 

dt 


MSJ) 


psi 

l+kl 


{d + y+£)I 


al 

co + f 


( 7 ) 


where the real, positive parameter a represents the maximal medical resources per unit time, 
and the real, positive parameter CO is the half-saturation constant. For simplicity, let the func¬ 
tions on the right-hand side of the equations in Q be /i and / 2 , respectively. Then, the 
equilibrium solutions of system (|^ are obtained by solving the following algebraic equa¬ 
tions: f\{S, I) = 0 and MS^ 7) = 0, from which the disease-free equilibrium can be easily 
obtained as Eq = {A/d, 0). For the infected equilibrium E = (S, /), S is solved from /i = 0 as 

5(7) = . Then, substituting S = S{I) into /2 = 0 yields a quadratic equation of 


the form 


{dk + P)I + d 


^{i) = £/f + dgi+^ = o. 


( 8 ) 


which in turn gives two roots: /i ,2 = where, ^ = {d + y + £){dk + p), 

SS = \^{dk -y P')CO d'^i^d y ~\~ £') ~y {^dk -f- /3)tt — /3A, ^ = ^(^d y £')(0 (X^d — PAco for 
system Q. Since all parameters take positive values, we have i?/ > 0. To get the two positive 
roots essential for backward bifurcation, it is required that <0 and ^ > 0. Noticing that 
/3, A, (0 > 0, we can see that the infection force, /3, the constant influx of the susceptibles, 
A, and the effect of medical treatment are indispensible terms for backward bifurcation. 
The number of positive infected equilibrium solutions changes from two to one when the 
value of C passes from negative to positive, which gives a critical point at C = 0, that is, 
[{d + y+£)(0 + a\d = pAco, which is equivalent to Rq = 

On the other hand, we may infer the emergence of backward bifurcation without solving 
the equilibrium conditions. If we do not consider the medical treatment term and remove 
it from system Q, that leads to system Q, which is a typical example of an SIR model studied 
by §. By setting the incidence function as /^{S, I) = we have /3(0, 7) = fsiS, 0) = 0; 
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Figure 1: Graphs of the ineidenee funetion /3 in system Q, (|^ and funetion /4 in system Q 
with respeet to 7, for whieh 5 = 50 has been used. The parameter values are ehosen as j 8 = 0.01, 
k = 0.01, a = 6, (0 = 1, d = 0.1, 7 = 0.01, e = 0.02, aeeording to P9l . The solid lines denote 
/s in (a) and /4 in (b), while the dashed ray lines in both graphs represent gi (7) = {d + 7 + e)7. 
(a) the ineidenee function / 3 (S, 7 ) = showing one intersection point with gi; and (b) the 
function / 4 ( 5 , 7 ) = showing two intersection points with gi. 


dfijSd) _ pi 


> 0 and 


PS 


> 0 for all 5, 7 > 0; and 


d^MSd) 


-IpkSil + 


~dr~ ~ T+ki ^ ^ 57 (i+/t /)2 ^ ^ ^ ^ w 

kl)~^ < 0 for all 5, 7 > 0. Therefore, the incidence function 7), satisfies the conditions 
given in Q. In particular, the function is concave, and can only have one intersection point 
with the line {d + 7 + e)I in the IS plane, as shown in Figure [^a). Thus, the uniqueness of 
the positive infected equilibrium implies that backward bifurcation cannot occur in this case. 
Moreover, according to the result in [j28]|, the uninfected and infected equilibria are globally 
asymptotically stable for Rq = < 1 and 7?o > 1, respectively. No complex dynamical 

behavior happens in system Q. 


In contrast, when we introduce the loss of the infectives due to medical treatment, the 
dynamics of system Q differ greatly from system Q. In particular, backward bifurcation 
emerges and complex dynamical behaviors may occur. To clarify this effect, we denote the 
function induced by 7 = 0 from as / 4 (S, 7) = Note that / 4 ( 5 , 7) is not an 

incidence rate. But, if we fix 5 = 5 > 0, there exist 0 < 7i < 72 < + 0 °, such that = 

- am(l+k7)2] > 0, V7 e (0, 72 ); and =-2k/35(l+k7)-3 + 

2am(m + 7)-3>0, V7e (0,7i), = 0, for 7 = 7i, <0, V7e ( 7 i, 72 ). Thus, 

/ 4 (S, 7) actually has a convex-concave ‘S’ shape, and may have two positive intersection points 
with the ray line, gi (7) = {d-\- Y+ e)7, in the first quadrant; see Figure[^b). These intersections 
contribute the two positive equilibrium solutions that are a necessary feature of backward 
bifurcation. 
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In summary we may conclude that the necessary terms which should be contained in sys¬ 
tem Q in order to have backward bifurcation are the constant influx A, the infection force j 8 , 
and the saturating medical treatment 


4.2 Backward bifurcation in the infection model with convex incidence 


Now we consider the 2-dimensional infection model Q which exhibits viral blips, studied 
in MM- The motivation for this model was a series of clinical discoveries indicating 
that viral infection can increase the density of a harmful chemical substance ||T9j|30l|39l|2^, 
thereby amplifying an associated biochemical reaction iHTI . and thus accelerating the infec¬ 
tion rate d. This cooperative phenomenon in viral infection is expressed by an increasing, 
saturating infectivity function: {B + According to the principle of mass action, the in¬ 

cidence function is then denoted as {B -f ^^)XY, which is a convex function with respect to 
the infectives’ density Y. 

To analyze the occurrence of possible backward bifurcation, we first examine the two equi¬ 
librium solutions from the following equations: 

AY AY 

fs (X, Y) = l-DX-iB+ = 0, fe{X, Y) = {B + = 0, (9) 

where all parameters A, B, C and D are positive constants. It is easy to find the uninfected equi¬ 
librium Eo = (Xo, ?o) = ( 5 , 0), whose characteristic polynomial has two roots: Ai = —D < 0, 
and A 2 = § — 1, which gives Rq = Consequently, Eq is stable (unstable) for < 1 (> 1). To 
find the infected equilibrium solution, setting f(,{X,Y) =0 yields Xi(y) = which 

is then substituted into / 5 (X, 7) = 0 to give the following quadratic equation: 

^ 5 ( 7 ) = (a+B)Y^ + {BC + D-A-B)Y + C{D-B)=0. (10) 


In order to have two real, positive roots, two conditions must be satisfied, that is, BC + D — 
A — B <0 and D — 5 > 0, or in compact form, 0 < D — B < A — BC. The condition D — B > 0 
is equivalent to 0 < = § < 1, which is a necessary condition for backward bifurcation. 

Moreover, the positive influx constant, having been scaled to 1, is a necessary term for the 
positive equilibrium of 7. Therefore, the positive influx rate term and the increasing and 
saturating infectivity function are necessary for backward bifurcation. 

In the rest of the subsection, we further examine the incidence function, 

AY 

MX,Y)^(B + —)XY, (11) 


without solving the equilibrium solutions. The incidence function // obviously satisfies the 
condition (j^l, as well as the condition (3b) since ■^f'j{X,Y) = [5-l-A7(74-C)“^]7 >0 and 
^fj{X, 7) = ACXY{Y + C)-^ + [B + A7(7 + C)-i]A > 0 for all X, 7 > 0. However, the 

second partial derivative of fiiX, 7) with respect to 7, -^fiiX, 7) = 2AC^X{X - 1 -C)^^ > 
0 for all X, 7 >0, showing that fqiX.Y) is a convex function with respect to the variable 
7. Consequently, / 7 (X, 7) can only have one intersection with g 2 (k) = 7, implying that 
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Figure 2: Graphs of the incidence functions /^{X.Y) and / 7 (F) for the parameter values A = 
0.364, B = 0.03, C = 0.823, and D = 0.057. The incidence functions are denoted by the 
solid lines, while the ray lines, determined by g 2 (F) = Y, are denoted by dotted lines: (a) the 
incidence function fqiX.Y), showing one intersection point with g 2 with an inset, with a fixed 
value X = 12.54; and (b) the incidence function fi{Y), showing two intersection points with 
an inset. 

only one equilibrium solution would exist if we only consider the second equation in Q, as 
shown Figure (a). However, when considering both conditions given in Q for equilibrium 
solutions, we will have two intersection points between // and gj. According to the first 
equation in (|^, that is / 5 (X, 7) = 0, we can use Y to express X in the equilibrium state as 
l(y) = (y + C)[(A+5)y2 + (5c+D)y+DC]-i. Substituting A(y) into//(X, y) in(|^, we 
obtain 

fn{Y)=Y[{A + B)Y+ BC][{A + B)Y^ + {BC + D)Y+ CD]-\ (12) 

and ^/v (y) = D[(A + 5)y2 + 2(A + B)CY + BC^] [{A + B)Y^ + {BC + D)Y + CD] ^>0 for 
all A, y > 0. However, the sign of = -2D[{A + BfY^ + 3C(A + BfY^ + 3(A + 

B)BC^Y + {if'C — AD)C'^] [{AAB)Y^ + {BC + D)Y +CD]^^, could alter at the inflection point 
from positive to negative as Y increases. Therefore, with appropriate parameter values, //(T) 
can have a convex-concave ‘S’ shape, yielding two intersection points with the ray line, g 2 {y), 
in the first quadrant of the X-Y plane, as shown in Figure |^(b). The above discussion, as 
illustrated in Figure implies that system Q can have two positive equilibrium solutions 
when < 1, and thus backward bifurcation may occur. 

Remark 1 Summarizing the discussions and results given in this section indicates that a dis¬ 
ease model with a convex-concave incidence function may lead to backward bifurcation, which 
in turn implies: (a) the system has at least two equilibrium solutions, and the two equilibrium 
solutions intersect at a transcritical bifurcation point; and (b) at least one of the equilibrium 
solutions is determined by a nonlinear equation. 
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5 Hopf bifurcation 


In the previous seetion, we studied baekward bifureation and established the neeessary eon- 
ditions for the oeeurrence of backward bifurcation in two models. In this section, we turn to 
Hopf bifurcation, since it typically underlies the change of stability in the upper branch of the 
infected equilibrium, the key condition in determining whether a model can exhibit oscillation 
or even recurrence. Again, we will present detailed studies for the two models. 


5.1 Hopf bifurcation in the infection model with concave incidence 


In this subsection, we study two cases of an infection model with concave incidence: sys¬ 
tem ([^ and Q. First, we discuss the equilibrium solutions and their stability by using the 
Jacobian matrix, denoted by 7, and examining the corresponding characteristic polynomial, 

F|y(L) =L2 + Tr(7)L-fDet(7). (13) 

Bifurcation analysis is conducted by choosing A as the bifurcation parameter. 

First, we consider the case without saturating medical treatment, system (|^. This system 
satisfies the three conditions in (3), and consequently, its uninfected equilibrium Eq = (^, 0) 

is globally asymptotically stable if Rq = < k while the infected equilibrium Ei = 

(rf^+j6Xrf+^+e) ) snisrges and is globally asymptotically stable if > 1- Therefore, 
for this case the system has only one transcritical bifurcation point at = 1 ^nd no complex 
dynamics can occur. 

Next, with the saturating treatment term, system (|7]) violates the conditions established for 
model Q, but leads to the possibility of complex dynamical behaviors. In fact, evaluating the 
Jacobian matrix Ji = 7| Q (Eq) at the uninfected equilibrium, Eq = (^, 0), yields the character¬ 
istic polynomial in the form of (13), denoted by P\j^ (L), with Tr(7i) = ^ , 

andDet(7i) = [—l5A + d^ + d£ + ^) = Tr{Ji)d — d^. This indicates that 
Det(7i) < 0 when Tr(7i) = 0, and thus Hopf bifurcation cannot occur from Eq. On the other 
hand, a static bifurcation can occur when Det(7i) = 0, that is, As = ^{d^+ d£ + where the 
subscript ‘S’ refers to static bifurcation. Therefore, Eq is stable (unstable) for A < As (> ^ 5 ), 

otRq < 1(> 1), with = liAd^^{d + Y+£ + [Bll. 

We will show that complex dynamical behaviors can emerge in system Q from the in¬ 
fected equilibrium Ei = (S, 7), where 7 is determined from the equation ^(7) = 0 in ([^. In 
the A-7 plane, the bifurcation diagram as shown in Eigure|^(l)-(4), indicates a turning point 
on the curve with appropriate parameter values, determined by both the quadratic equation Q 
and the relation ^ = 0, which is equivalent to ^ = 0. Solving ^ = 0 yields 

the turning point of 7, denoted by 7^ (‘T’ means turning), taking the form 


h = 7: 


PAt 


(0 


a 


2 (7k-f"/3) (7-(-c) dk-\~P d-\-£ 

where Aj is obtained from = 0, see ([^. Thus, when 7^ > 0 (< 0), the turning point of 

the quadratic curve appears above (below) the 7-axis, meaning that backward bifurcation oc¬ 
curs for 7 > 0 (< 0). Evaluating the Jacobian matrix at the infected equilibrium Ei, and further 
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denoting it as 72 = ■^|(( 7 l)(Ei), we obtain the characteristic polynomial in the form of (13 
Tr(72) = ai 1 /[((0 + lY{kl + 1 ){dkl + /3/ + 7)] and Det(72) =a 2 i/[{ 0 ) + lY{kl+l){dkl- 


, with 

/3/+ 


d)], whereaii =aia — a\b anda 2 i =a2a — a2b, withai^ = PA{(0 + iy anda2b = da\b, andai^ 
and a 2 a only contain positive terms (their expressions are omitted here for brevity). Therefore, 
we can rewrite Det( 72 ) = ^d/[{co + I)^{kl + l){dkl + + d)]. Determining whetheraHopf 
bifurcation can occur from E is equivalent to finding whether Det(72) remains positive when 
Tr(72) = 0. Ignoring the positive factors in the following subtraction yields 


, ,,, Tr(72)-Det(72)/J 1 1 

{co + lY{kI+l){dkI + pi + d) d""^^ 

where/zi(/) = ^{dkl + + d)[{kl +l)d^{a}+I)^ — P£l{(0 + I)^ — OcPcoI]. Thus, whenaia = 

0, ^a 2 a and hi (I) have opposite signs, implying that when Tr(72) = 0, Det(72) could be positive 
only if hi{I) is negative. Therefore, the necessary condition for system Q to have a Hopf 
bifurcation from the infected equilibrium Ei is that hi (I) is negative. 


In the remaining part of this subsection, we demonstrate various dynamics which may 
happen in system (|7]) with different parameter values of k, as shown in Table Taking other 
parameter values as a = 6 , CO = 7, e = 0.02, 7 = 0.01, /3 = 0.01, and d = 0.1, and solving 
the two equations Tr( 72 ) = 0 and ^{I) = Qm^ gives the Hopf bifurcation point candidates, 
(A//, Ih), for which hi (7//) < 0. Since the formula for the transcritical bifurcation point A 5 
has no relation with k, (A 5 , Is) = (9.87, 0) is a fixed value pair in Table Bifurcation dia¬ 
grams and associated numerical simulations are shown in Figure corresponding to the five 
cases given in Table [I] The blue lines and red curves represent the uninfected equilibrium Eq 
and infected equilibrium Ei, respectively. The stable and unstable equilibrium solutions are 
shown by solid and dashed lines/curves, respectively. Backward bifurcation occurs in Cases 
1, 2, and 3 (see Table [^, which are illustrated by the corresponding bifurcation diagrams in 
Figuresj^l), (2), and (3), respectively. For Cases 1 and 2, only one Hopf bifurcation occurs on 
the upper branch of the infected equilibrium Ei, and this bifurcation point exists at the critical 
point Ah < A 5 for Case 1 and Ah > As, for Case 2. For Case 1 with A = 9.78, the simulated 
time history converges to Eq with initial condition IC= [93.6, 0.44], shown in Figure |^la), 
but converges to Ei with initial condition IC= [46.8, 10], shown in Figure|^lb). This clearly 
indicates the bistable behavior when Ah < As, and an overlapping stable region for both Eq 
and El exists (see Figure [^1)). The recurrent behavior for Case 2 is simulated at A = 9.87 
with IC= [50, 5], shown in Figure |^2a). For Case 2, Ah > As, and an overlapping unstable 
parameter region for both Eq and Ei occurs between As and Ah (see Figure[^2)). For Case 3, 
two Hopf bifurcations occur on the left side of As, and a stable part in the upper branch of Ei 
exists when A passes through the critical value A = A 5 . In this case, although backward bifur¬ 
cation still exists and the turning point is also located above the A-axis, giving two branches 
of biologically feasible Ei, only regular oscillating behavior is observed. The simulated time 
history is conducted at A = 10, with initial condition IC= [50, 2], shown in Figure |^3a). For 
Case 4, only forward bifurcation occurs in the biologically feasible region, and the turning 
point for backward bifurcation moves down to the fourth quadrant, that is, negative backward 
bifurcation occurs in this case. The whole upper branch of Ei in the first quadrant is stable, 
therefore, no oscillations (or recurrence) can happen. Finally, further increases to the value of 
k change the shape of the red curves, as shown in Figure |^5), which again indicates that no bi¬ 
ologically meaningful backward bifurcation or oscillations can occur. Note that in Figure [^5) 
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Table 1: Dynamics of system Q for different values of k, with a = 6, CO = 7, £ = 0.02, 
7 = 0.01, j 8 = 0.01, d = 0.1, and a fixed transcritical bifurcation point (A 5 , 75 ) = (9.87, 0). 


Case 

k 

(Ar, h) 

hi{I)<0 

(A/f, Ih) 

Dynamics 

Notes 

1 

0.001 

(9.48,4.57) 

7 G [ 1 . 72 , 00 ] 

(9.73, 10.28) 

Bistability 

A// < A 5 

2 

0.01 

(9.71,2.82) 

7 G [1.76, 00 ] 

(9.96, 8.00) 

Recurrence 

A/f > A 5 

3 

0.02 

(9.85,0.84) 

7 G [1.82, 00 ] 

(9.88, 2.09), 
(10.14,5.62) 

Oscillation 

Two Hopf 
critical points 

4 

0.027 

(9.86, -0.65) 

7 G [1.85,30.65] 

No Hopf 

No 

oscillation 

Negative 

backward 

bifurcation 

5 

0.05 

No Turning 

7 G [2.01,15.03] 

(6.18,-22.15) 

No 

oscillation 

No backward 
bifurcation 


a Hopf bifurcation point exists on the lower branch of the equilibrium solution, which is bi¬ 
ologically unfeasible since it is entirely below the horizontal axis. In conclusion, interesting 
dynamical behaviors can emerge in system Q if backward bifurcation occurs. 


5.2 Hopf bifurcation in the infection model with convex incidence 


In this subsection, we return to system Q, that is, the 2-dimensional HIV model with convex 
incidence derived in [|47ll4^ . and analyze the various dynamical phenomena which system Q 
could possibly exhibit. To achieve this, we set B as the bifurcation parameter, and A as a control 
parameter; the bifurcation analysis will be carried out for various values of A. Also, simulated 
time histories are provided to illustrate the dynamical behavior predicted in the analysis. 

We first consider the uninfected equilibrium Eq = which has two eigenvalues. One 

of them, given by = —D, is always negative. The other one is ^ — 1. Thus, 

depending upon the relation between B and D, X 2 \e^ = 0 gives a static bifurcation at Bs = D 
(or 7?o = § = 1), which is further proved to be a transcritical bifurcation. Here the ‘S’ in 
subscript stands for static bifurcation. Therefore, Eq is stable when B < D {or Rq < 1), loses 
its stability and becomes unstable when B increases to pass through Bs = D, that isB > D (or 
^0 > 1), and no other bifurcations can happen. 

Next, we examine the infected equilibrium Ei = (X, 7). Since A(y) = ^ 

determined by the quadratic equation ([T^, which gives the turning point (5^, Yj) as 


Bj 


-A + D + 2VACD A + B-BC-D 

C+1 ’ “ A+B 


where ‘T’ in the subscript stands for turning bifurcation. 


We perform a further bifurcation 
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Figure 3: Bifurcation diagrams and simulations associated with the five cases given in Table 
demonstrating various dynamical behaviors. 
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analysis on its corresponding characteristic polynomial ( fTS] ), which takes the form 

P\p iX.Y) = + Y , - ^ — — — T7 - + 77-7—~—77-7? whcrc 

^ [(A+5)y+5c](y+c) [(A+5)y+5c](y+c)’ 

= {A + BYY^ + {2BC + D){A + B)Y^ + {B^C'^+ACD + 2BCD-AC)Y + BC^D, 
aia = (A + 5)2y3 + 2{A+ B)BCY^ + {B^C-AD)CY. 


Therefore, the sign of the subtraction between the trace and determinant is determined by 
h2{Y) = a\a — a2a = D{A + B)Y^ + [ 2 CD{A + B) — AC]y +BC^D. Here the equilibrium solu¬ 
tion of y and other parameters satisfy the quadratic equation ( [T^ , which leads to an explicit 

expression, given by 5 = ~ y2+^(c-'i)y-c^ ' Substituting B = B into h2{Y), we obtain 

, ,,,,, [AC(D-\)-D^]Y^-[AC(D-\) + 2CD'^]Y-C^D'^ 

h2{Y)\B=B = aia-a2a= - ^7 -^-■ 


Hopf bifurcation may occur when the trace is zero, while the determinant is still positive. This 
implies h 2 {Y) < 0, which is possible with appropriately chosen parameter values. Hence, by 
solving aia = 0 in ( fH] ) together with the quadratic equation ( [TO] ), we get two pairs of points 
denoted by (5/,i, Y^i) and (5 /j2 , y/, 2 ), which are candidates for Hopf bifurcation. Then vali¬ 
dating the above two points by substituting them back into the characteristic polynomial ( [T4| ), 
respectively, we denote the Hopf bifurcation point as (Bh^Yh) if this validation confirms their 
existence. According to (Yu et ah, submitted for publication), Hopf bifurcation can happen 
only from the upper branch of the infected equilibrium Ei. 


The various dynamical behaviors which may appear in system Q have been classified in 
Tablej^for different values of the parameter A, with fixed values of C = 0.823 and D = 0.057. 
Thus, the transcritical bifurcation point is fixed for all cases: Bs = D = 0.057 and y 5 = 0. The 
two solutions Bpi and 5 /j2 are solved from the two equations (14) P\b^{X,Y) = 0 and (10) 
^ 5 {Y) = 0, respectively. They become a Hopf bifurcation point only if their corresponding Y 
values (y /,1 and Yh 2 , respectively) are in the range such that h 2 {Y) < 0. Otherwise, system (|^ 
has a pair of real eigenvalues with opposite signs at (5/,i, Yhi) or (5/,2, Yh 2 ), which is denoted 
by the superscript (which is actually a saddle point) in Table while the Hopf bifurcation 
point is denoted by the superscript ‘77’ in Table 


Next, we further examine the direction of the Hopf bifurcation, that is, check whether it is a 
supercritical or subcritical Hopf bifurcation. Since the Jacobian matrix of the system evaluated 
at the Hopf bifurcation point has a pair of purely imaginary eigenvalues, the linearized system 
(|^ does not determine the nonlinear behavior of the system. Therefore, we take advantage 
of normal form theory to study the existence of the limit cycles bifurcating from the Hopf 
bifurcation point as well as their stability. As mentioned earlier, Hopf bifurcation can only 
occur from the upper branch of the infected equilibrium E\, therefore we first transform the 
fixed point Ei to the origin by a shifting transformation, and, in addition, make the parameter 
transformation B = Bh + /i; the Hopf bifurcation point is thus defined as /i = /i// = 0. Then the 
normal form of system (Q near the critical point, /i = /i// = 0, takes the form up to third-order 
approximation: 


r = d lJ.r +ar^ + 9 = G)c + cjl-\-br^ + (15) 
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Table 2: Parameter values taken to illustrate various dynamics of system Q. The fixed 


transcritical bifurcation point: {Bs, Ys 

= (0.057, 0) 

Case 

A 

(5r, Yt) 

/?2(t) <0,ye 

Yhi) 

1 

0.80 

(-0.1950, 0.5850) 

(0.0036, 0.9830) 

(0.0355, 0.8725)" 

2 

0.71 

(-0.1580, 0.5660) 

(0.0040, 0.9800) 

(0.0539, 0.0038)* 

3 

0.60 

(-0.1140, 0.5380) 

(0.0048, 0.9769) 

(0.0540, 0.0045)* 

4 

0.07 

(0.0557, 0.0909) 

(0.0476, 0.8030) 

(0.0560, 0.0470)* 

5 

0.06 

(0.056558,0.05581) 

(0.0574, 0.7700) 

(0.056559, 0.0574)" 

6 

0.05 

(0.05697, 0.01442) 

(0.0724, 0.7232) 

(0.0574, 0.0741)" 

7 

0.04 

(0.0569, -0.0358) 

(0.0986, 0.6507) 

(0.0592, 0.1071)" 

8 

0.03 

(0.0559, -0.0994) 

(0.1611,0.5149) 

— 

Case 

A 

{Bhi, Yhi) 

Dynamics 

Notes 

1 

0.80 

(0.054, 0.0034)* 

Unstable limit cycle, 
Bistable 

Bhi < Bs 

2 

0.71 

(0.0574, 0.8650)" 

Recurrence 

Bhi > Bs 

3 

0.60 

(0.0819,0.8530)" 

Recurrence 

Bhi > Bs 

4 

0.07 

(0.1015,0.5612)" 

Recurrence 

Bhi > Bs 

5 

0.06 

(0.0961,0.5225)" 

Recurrence 

Bhi <Bs < Bhi 

6 

0.05 

(0.0894, 0.4701)" 

Recurrence 

Bhi <Bs < Bhi 

7 

0.04 

(0.0806, 0.3897)" 

Oscillation 

Bhi < Bs < Bhi, 
Yt<0 

8 

0.03 

— 

Ei stable 

Yr<0 
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Table 3: Classification of Hopf bifurcations based on the normal form ([T5|). 


Class 

Stability of r = 0 

Stability of r^ = -^ 

Hopf bifurcation 

11 <0 

li>0 

11 <0 

11 >0 

(a): d > 0, a> 0 

stable 

unstable 

unstable 

- 

subcritical 

(b): d > 0, a <0 

stable 

unstable 

- 

stable 

supercritical 

(c): d <0, a> 0 

unstable 

stable 

- 

unstable 

subcritical 

(d): d < 0, a < 0 

unstable 

stable 

stable 

- 

supercritical 


Table 4: Classification of Hopf bifurcations appearing in Table 


Case 

A 

Hopf bifurcation 

point {Bh, Yh) 

d 

a 

Stability of 
limit cycles 

Tabled 

class 

1 

0.8 

(0.0355, 0.8725) 

-1.0722 

0.2114 X 10~^ 

Unstable 

(c) 

2 

0.71 

(0.0574, 0.8650) 

-1.0726 

0.1424 X 10-^ 

Unstable 

(c) 

3 

0.6 

(0.0819,0.8530) 

-1.0733 

0.6755 X 10-^ 

Unstable 

(C) 

4 

0.07 

(0.1015,0.5612) 

-1.0307 

-0.8791 X 10^^ 

stable 

(d) 

5 

0.06 

(0.056559, 0.0574) 

884.27 

-0.1019 

Stable 

(b) 



(0.0961,0.5225) 

-1.0079 

-0.8613 X 10-^ 

Stable 

(d) 

6 

0.05 

(0.0574, 0.0741) 

18.232 

-0.3145 X 10-^ 

Stable 

(b) 



(0.0894, 0.4701) 

-0.9629 

-0.8457 X 10-^ 

Stable 

(d) 

7 

0.04 

(0.0592,0.1071) 

4.7242 

-0.1577 X 10^^ 

Stable 

(b) 



(0.0805, 0.3897) 

-0.8437 

-0.8438 X 10-^ 

Stable 

(d) 


where r and 0 represent the amplitude and phase of the motion, respectively. The first equation 
of ( 15 1 can be used for bifurcation and stability analysis, while the second equation of ( [T5| ) can 
be used to determine the frequency of the bifurcating periodic motions. The positive (Oc in the 
second equation of ( [T5| ) is the imaginary part of the eigenvalues at the Hopf bifurcation point. 
The parameters d and c can be easily obtained from a linear analysis, while a and b must be 
derived using a nonlinear analysis, with the Maple program available in, say, fl46ll . 

Note that the infected equilibrium Ei is represented by the fixed point r = 0 of system 
(151, while the nonzero fixed point r > 0 (satisfying is an approximate solution for 

a limit cycle or periodic orbit. The periodic orbit is asymptotically stable (unstable) if a < 0 
{a > 0), and the corresponding Hopf bifurcation is called supercritical (subcritical). According 
to the Poincare-Andronov Hopf Bifurcation theorem [l44ll . for /i sufficiently small, there are 
four possibilities for the existence of periodic orbits and their stability, which are classified 
in Table based on the four sets of the parameter values in the normal form ( [T5] ). Then we 
use the results presented in Table with a nonlinear analysis based on normal form theory to 
classify the Hopf bifurcations appearing in Table and the results are shown in Table 

To illustrate the analytical results given in Tables and we provide the bifurcation di¬ 
agrams in Figures 1^ (l)-(8). These figures depict the uninfected equilibrium Eq and the in¬ 
fected equilibrium Ei in blue and red, respectively. The solid and dashed lines differentiate 
stable and unstable states of the equilibrium solutions. The bifurcation points on the equi¬ 
librium solutions are highlighted by solid black dots. Moreover, ‘Transcritical’, ‘Turning’, 
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‘Hopfsub’j and ‘Hopfguper’, are used to denote Transcritical bifurcation. Turning point, sub- 
critical Hopf bifurcation, and supercritical Hopf bifurcation, respeetively. Simulated time 
histories are used to validate the analytieal results, and to show different dynamieal behav¬ 
iors in each case listed in Tables and Subcritical Hopf bifurcation occurs in Cases 
1-3, shown in Figures (l)-(3). A = 0.8 is used in Figure (1) for Case 1. Choosing 
B = 0.036, we have Eq = [17.1282566, 0.023689] and Ei = [2.233533, 0.8726886]. The sim¬ 
ulated solution converges to Eq or Ei, with initial condition taken as IC^/ = [17.13,0.024] 
or ICc = [2.233, 0.873], shown in Eigures (Id) and (Ic), respectively. Eigures (la) and 
(lb), on the other hand, show the unstable limit cycle bifurcating from the subcritical Hopf 
bifurcation with ICc = [2.233, 0.873]. 

Eigure|^(2) corresponds to Case 2 with A = 0.71. Choosing B = 0.0572 G [^5, Bh] yields 
recurrence, independent of the initial conditions, see, for example, the result given in Eigure]^ 
(2b) with ICi, = [2.4, 0.5]. However, for B = 0.06 > Bh, the simulated time history converges 
to El, with an initial condition close to Ei, such as ICq = [2.4, 0.6] as shown in Eigure|^(2a); 
or shows recurrence with an initial condition far away from Ei, such as ICc = [2-4, 0.4], as 
shown in Eigure|^(2c). 

Eigure|^(3) plots the result for Case 3 with A = 0.6, and shows a broader region between 
the transcritical and Hopf bifurcation points, associated with a larger recurrent region. Re¬ 
currence occurs independent of the initial conditions for B = 0.083 G [fi^, Bh], giving Eq = 
[12.048, 0] and El = [2.576, 0.852], as shown in Eigures|^(3a) and (3b), withICa = [2.7, 0.84] 
and IC;, = [14, 0.1], respectively. But if we choose B = 0.07 > Bh, we have Eq = [14.286,0] 
and El = [2.67, 0.8478]. The time history converges to Ei with ICc = [2.6, 0.8], or shows 
recurrence with IC^ = [2.6, O.I], as shown in Eigure|^(3c) and (3d), respectively. 

Supercritical Hopf bifurcations occur in Cases 4-7, as shown in Eigures |5(4)-(7). Eig- 
ure|^(4) depicts the result for Case 4 with A = 0.07. Only one supercritical Hopf bifurcation 
happens in this case, and gives a large recurrent parameter region between the transcritical 
and Hopf bifurcation points. Although the simulated recurrent behavior does not depend on 
initial conditions, the recurrent pattern will fade out with the growth of the value of B from 
the transcritical point to the Hopf bifurcation point, see Eigures]^ (4a) and (4b) with the same 
^Ca^b = [8, 0.1], but different values of B: B = 0.06 and B = 0.09, respectively. 

Eigure|^(5) shows the result for Case 5 with A = 0.06. A transcritical bifurcation hap¬ 
pens between two supercritical Hopf bifurcations. The recurrent region still starts from the 
transcritical point and independent of the initial conditions, but is narrower than that shown in 
Eigure|^(4). The simulated recurrent behavior for this case is conducted at IC= [12, 0.1] and 
B = 0.06. Eigure|^( 6 ) corresponds to Case 6 with A = 0.05, and two supercritical Hopf bifur¬ 
cations occur on the right side of the transcritical bifurcation point, which makes the recurrent 
region even narrower and the recurrent pattern less obvious, as shown in the simulated time 
history with IC= [10, 0.1] and B = 0.06. Negative backward bifurcations occur in Cases 7 and 
8 , as shown in Eigure|^(7) and ( 8 ). Although two Hopf bifurcations are still present in Case 
7, see Eigure|^(7), only a regular oscillating pattern exists. Eor Case 8 , no Hopf bifurcation 
happens in the biologically feasible part of Ei, and therefore no more interesting dynamics 
occur. 

In general, backward bifurcation, which occurs above the horizontal axis, is much more 
likely to induce Hopf bifurcation. A Hopf bifurcation can only occur along the upper branch 
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of since Eq only changes its stability at a transcritical bifurcation point, and any point 
on the lower branch of E\ is a saddle node (Yu et ah, submitted for publication). Moreover, 
Hopf bifurcation can lead to a change in the stability of the upper branch of the infected 
equilibrium E\. Thus the system further develops bistable, recurrent, or regular oscillating 
behavior, corresponding to Cases 1 — 7 in Tables|^and|^ and in Figures|^(l)-(7). In particular, 
bistability happens when both equilibria Eq and E\ share a stable parameter region, see Case 
1 in Tableand Figure]^ (1). 

As for recurrent behavior, we observe that recurrence is more likely to happen if the fol¬ 
lowing two conditions are satisfied for the upper branch of Ei: (1) the equilibrium remains 
unstable as the bifurcation parameter increases and crosses the trancritical point, where Eq and 
El intersect, such that the two equilibria share an unstable parameter range; and (2) at least one 
Hopf bifurcation occurs from Ei. As shown in Cases 2-6 in Table and the corresponding 
Figures|^(2)-(5), the common recurrent parameter region for both subcritical and supercritical 
Hopf bifurcations starts beside the transcritical point, and is located entirely in the unstable 
parameter region of Eq and Ei. The simulated recurrent pattern becomes more pronounced 
if the value of the bifurcation parameter is close to the transcritical point, but approaches an 
oscillatory pattern as the parameter diverges from the transcritical point, as shown in Figure]^ 
(4a) and (4b). In this common recurrent parameter region, recurrence occurs independent of 
initial conditions; see Figures (3a) and (3b). In addition to the common recurrent region, 
for subcritical bifurcation, seen in Table for Cases (2) and (3) and Figures (2) and (3), 
recurrence may also appear on the stable side of the subcritical Hopf bifurcation point with an 
initial condition close to Ei. Moreover, the subcritical Hopf bifurcation and the transcritical 
point should be close to each other for a clear recurrent pattern. When this is not the case, 
the periodic solutions show a more regular oscillating pattern, as compared in Figures (2c) 
and (3d). Although two Hopf bifurcation points occur in Tablej^for Case 5, see Figure|^(5), 
the transcritical point is located inside the unstable range of the upper branch of £ 1 , between 
the two Hopf bifurcation points. A recurrent pattern still characterizes the dynamical behavior 
in this case. However, if the unstable range of Ei, between the two Hopf bifurcation points, 
is located entirely in the unstable range of Eq, and moves further away from the transcritical 
point, the recurrent motion gradually becomes a regular oscillation, as shown in Figures (6) 
and (7). 

Summarizing the results and discussions presented in the previous two sections, we have 
the following observations. 

1. Due to the fact that Eq only changes its stability at the transcritical bifurcation point, and 
the fact that any point on the lower branch of Ei is a saddle node, Hopf bifurcation can 
only occur from the upper branch of Ei. A Hopf bifurcation may result in convergent, 
recurrent, bistable, or regular oscillating behaviors. 

2. Backward bifurcation gives rise to two branches in the infected equilibrium Ei. Hopf 
bifurcation is more likely to happen when the turning point of the backward bifurcation 
is located on the positive part of the equilibrium solution in the bifurcation diagram, as 
shown in Figures |^(2)-(6). This means that we have two biologically feasible infected 
equilibria, which is essential to observe bistability, as shown in Figure|^(l). 

3. However, if the turning point on the infected equilibrium Ei, or the backward bifurcation 


18 


moves down to the negative part of a state variable in the bifureation diagram, that is, 
negative baekward bifureation oeeurs, then Hopf bifureation is very unlikely to happen. 
Although Figure (7) shows an exeeptional ease, the parameter range for sueh a Hopf 
bifurcation is very narrow. 

4. The bifurcation diagram for system Q with A = 0.03, shown in Figure|^(8), is a typical 
model with negative backward bifurcation. Such negative backward bifurcation may 
occur in higher-dimensional systems. However, by considering more state variables, 
which make the system more complicated, Hopf bifurcation can happen in the upper 
branch of the negative backward bifurcation. We will discuss this possibility in more 
detail in the next section by examining an autoimmune disease model. 

The results obtained in this section suggest the following summary. 

Remark 2 If a disease model contains a backward bifurcation on an equilibrium solution, 
then as the system parameters are varied, there may exist none, one or two Hopf bifurcations 
from the equilibrium solution, which may be supercritical or subcritical. If further this equi¬ 
librium has a transcritical bifurcation point at which it exchanges its stability with another 
equilibrium, then recurrence can occur between the transcritical and Hopf bifurcation points 
and near the transcritical point, where both equilibrium solutions are unstable, and bistability 
happens when Hopf bifurcation makes a shared stable parameter region for both equilibria. 


6 Negative backward bifurcation in an autoimmune disease 
model 


In the previous section, we examined three cases of negative backward bifurcation: Table 
Case 4 for system (|^ and Table|^Case (7) and (8) for system Q. The analytical and numerical 
results showed that solutions typically converge to the infected equilibrium in these cases, 
and the parameter range for Hopf bifurcation is very limited. As a result, negative backward 
bifurcation tends to give no interesting behavior. In this section, however, we shall explore 
an established autoimmune model [|T]| in which negative backward bifurcation occurs. We 
demonstrate that after modification, the autoimmune model can also exhibit recurrence. 

The autoimmune model [[B takes the form 


cL4 

dt 

dRn 


dt 

dE 

dt 

dG 

dt 


fvG- {oiRn + bi)A- PaA 
-- {TtxEA^)A-pnRn 
XeA — PeE 
jE-vG-PgG, 


(16) 


where mature pAPCs (A) undergo maturation by intaking self-antigen (G), at rate /v, and are 
suppressed by specific regulatory T cells, TReg cells at rate op b\ represents additional 
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( 1 ) 


( 2 ) 





Figure 4: Dynamical behaviors of system Q corresponding to eight cases listed in Tablej^and 
1^ All insets are simulated time histories of 7 vs. t. The yellow areas fading to white show 
regions in which recurrent behavior occurs and fades to regular oscillations. 
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Figure 4: Dynamical behaviors of system Q corresponding to eight cases listed in Tablej^and 
1^ All insets are simulated time histories of 7 vs. t. The yellow areas fading to white show 
regions in which recurrent behavior occurs and fades to regular oscillations. 
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non-specific background suppression. The TReg cells are activated by mature pAPCs at a rate 
proportional to the number of auto-reactive effector T cells {E) at rate TZ\, and by other sources 
at rate /3. Active auto-reactive effector T cells {E) come from the activation process initiated 
by mature pAPCs, at rate Xe, then attack healthy body tissue and release free self-antigen (G) 
at rate 7 , which is ready for mature pAPCs to engulf; the antigen engulfing rate is v. The death 
rates of the populations A, Rn, E, and G are denoted by /i„, /i£, and /ic, respectively. 

Following the steps described by Zhang et al. (submitted for publication), system ( fT^ , can 
be reduced via quasi-steady state analysis to a 2 -dimensional system: 


dA 


dt 

dRn 


dt 


fvy^E 

l^Eiy + jic) 

jJ-E 


-bi- jUa]A - OiRnA, 


■ l-^nEn 


(17) 


For simplicity, we set a = —b\—jlA and b = For the stability and bifurcation 

analysis, we choose Xe as the bifurcation parameter. System ( flT] ) has a disease-free equilib¬ 
rium Eq = (0, 0), which is stable if a > 0 or ; and unstable if a < 0 or < 

{b\EiiA){vAiLc)^E ^ ^ static bifurcation occurs on Eq when a = 0 or A^ = _ 

The disease equilibrium is given by £1 = (A, ^„), in which R,^ = and A is given by 

the roots of the following equation. 


/ 8 (A) = bOiA^ + j5aiA- iina. 


(18) 


Equation ( [18] ) has two roots with negative signs if a < 0, with opposite signs if a > 0, and 
only one zero root if a = 0. This means that a negative backward bifurcation is possible in 
system ( [T7] ) with proper parameter values. We further examine the characteristic equation 
at El, which shares the same form as equation (13), with Tr(y|£j) = -^{baiA^ + j5oiA + 


— ajin) ■= «ii and Det(7|£j) = 3baiA^ -I- 2/3aiA — ajin := < 212 . Solving / 8 (A) = 0 and 
ai 2 = Det(7|£j) = 0, gives the static bifurcation point of Ei at (A, a) = (0, 0) or (A, Xe) = 
(0, )^ which is a transcritical bifurcation point between Eq and Ei. Moreover, 

Hopf bifurcation can happen if and only if / 8 (A) = 0 and an = Tr(7|£) = 0, which can be 
satisfied only if jU„ = 0. This implies that the positive branch of £1 is stable for any positive 
values of /i„. Thus, this model cannot exhibit recurrence, bistability, or even regular oscillation. 
The same conclusion was obtained in Zhang et ah, (submitted for publication) for the original 
4-dimensional model ( [T^ . 

However, a recent experimental discovery [|3]| has revealed a new class of terminally dif¬ 
ferentiated TReg cells. As described in detail in Zhang et ah, (submitted for publication), 
introducing this cell population, denoted Rj, into the model yields the full system 


^ = /vG - CTi {Rn + dRd)A - {bi -f iia)A 
^ = {niE+^)A-^nRn-^Rn 
^=C^Rn-lldRd 
^ = A^A — iieE 
f = 7E-vG-ilGG 


22 

















and quasi-steady state analysis then yields a redueed 3-dimensional model in the form 


Again, here Ae is chosen as the bifurcation parameter for stability and bifurcation analysis. It is 
easy to show that system (|^ still has a disease-free equilibrium Eq as (A, R^, Rj) = (0, 0, 0), 
and a disease equilibrium E\ as (A, R^ Rd), where Rd = Rn = ^ 

determined from the following quadratic equation: 

MA) = n, XeA^ + /SjieA + [-fjvlE + (bi+HA)(tia + v)ji«] , 

which gives two negative roots if Ae < Aes = and two roots with opposite 

signs when Ae > Aes- The critical point is determined by Ae = Aes, which is actually the 
intersection point of Eq and Ei. The two equilibrium solutions exchange their stability at Aes, 
leading to a transcritical bifurcation at {A, Ae) = ( 0 , ^£ 5 ). Note that the negative backward 
bifurcation still happens in system Moreover, a Hopf bifurcation occurs from the upper 
branch of £ 1 , giving rise to oscillation and recurrence. 

Realistic parameter values have been obtained in Zhang et ah, (submitted for publication), 
and are given as follows: 

/=lxl0~4, v = 0.25x10-2, c7i =3x10-6, = 0.25, ^tA = 0.2, ;ri= 0.016, 

j 8 =200, jin = 0.1, jU£=0.2, 7=2000, I1g = 5, jid = 0.2, c = 8, d = 2, <^=0.025. 

For the above parameter values, the Hopf critical point is obtained at [Ah, Aeh) = (5.6739, 
1691.6414), while the turning point is at {Aj, Aet) — (—1.4205, 879.9848), and the trans¬ 
critical bifurcation point is at {As, Aes) = (0, 900.45). These three bifurcation points and the 
stability of equilibrium solutions are shown in the bifurcation diagram given in Figure [^a), 
and the simulated recurrent time history is plotted in Figurej^b) for A£ = Aeh + 1000. 

In summary, when negative backward bifurcation occurs, that is, the turning point is lo¬ 
cated in the negative state variable space, less complex dynamical behavior will be present. 
Hopf bifurcation in a biologically feasible area does not happen in the reduced 2-dimensional 
system nor in the original system ( [T^ (Zhang et ah, submitted for publication). How¬ 
ever, if we increase the dimension of the system, Hopf bifurcation and complex dynamical 
phenomena can emerge, as shown in our results for system Q. 


7 Conclusion 

In this paper, we first review previous work on a reduced 2-dimensional infection model with 
a concave incidence rate ||28l . The authors proved that the disease equilibrium will emerge 
and be globally stable when the basic reproduction number Rq is greater than 1. This means 
that no complex dynamical phenomenon can occur in such models. However, by adding an 
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(a) 


(b) 




Figure 5: Dynamics of system ( [T7| ): (a) bifurcation diagram; and (b) simulated time history 
for Xe = Xeh + 1000. 


extra saturating treatment term to this simple 2-dimensional infection model, the resulting 
system Q considered in P9l can exhibit backward bifurcation, which increases the parameter 
range for Hopf bifurcation, which in turn leads to recurrent, bistable and regular oscillating 
behaviors. 

Instead of adding an extra term, a 2-dimensional infection model with a convex incidence 
function can likewise show rich dynamics due to the occurrence of backward bifurcation, giv¬ 
ing rise to two types of Hopf bifurcation. Biologically, a convex incidence rate implies that 
existing infection makes the host more vulnerable to further infection, showing a cooperative 
effect in disease progression. From the view point of mathematics, the convex incidence func¬ 
tion enables backward bifurcation to occur on the positive branch of the disease equilibrium 
solution, which further generates Hopf bifurcation. The location and direction of Hopf bi- 
furcation(s), determined by parameter values, can further give rise to bistable, recurrent, and 
regular oscillating behaviors. 

Cooperative effects also occur during the progression of autoimmune disease. However, 
for an autoimmune model with negative backward bifurcation, in which the turning point 
is located on the negative state variable space, the biologically feasible parameter range in 
which Hopf bifurcation may occur is limited. By introducing an additional state variable to 
the autoimmune model, recurrent phenomenon are once again observed. 
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